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SECTION  1 
INTRODUCTION 


One  of  the  most  difficult  fluid  dynamic  analysis  problems,  in  the 
design  of  chemical  lasers,  is  the  computation  of  the  mixing  between 
the  primary  flow  and  various  injected  secondary  flows.  The  exact 
calculation  of  the  mixing  process  requires  the  solution  of  the 
three-dimensional  (3-D)  Navier-Stokes  equations  for  a  chemically 
reacting,  viscous,  compressible  flow  over  a  very  fine  mesh.  Due  to 
the  tight  mesh  requirements  and  the  need  to  solve  large  numbers  of 
species  equations,  exact  solutions  require  hundreds  of  hours  of 
computer  time  per  case,  even  on  CRAY  class  machines.  As  a 
consequence,  chemical  laser  mixing  analyses  have  been  done  almost 
exclusively  using  the  two-dimensional  (either  axisymmetric  or 
planar  2-D)  parabolic,  thin  shear  layer  approximation  (Ref.  1)  . 
When  the  secondary  flow  is  injected  parallel  to  the  primary  flow, 
the  2-D  approach  is  valid.  However,  in  most  instances,  the 
sidewall  injection  concept  is  used  wherein  the  secondary  flow  is 
injected  into  the  primary  at  some  large  angle.  In  this  instance, 
the  bending  of  the  injected  jet  in  the  direction  of  the  primary 
flow  induces  a  large  axial  vorticity  (i.e.,  flow  swirl)  component 
into  the  jet.  Since  the  2-D  equations  have  zero  axial  vorticity  by 
definition,  there  is  every  reason  to  believe  that  the  2-D  approach 
will  underpredict  the  rate  at  which  the  two  flows  will  mix. 

The  rate  at  which  mixing  takes  place  is  dependent  upon  the 
magnitude  of  the  convective  velocities  transverse  to  the  flow 
direction  and  the  magnitude  of  the  diffusional  velocities.  With 
the  2-D  shear  layer  equations,  the  transverse  convective  velocity 
is  small  and  the  choice  of  a  model  for  the  diffusional  velocities 
(i.e.,  binary  diffusion  versus  multicomponent  diffusion)  has  a 
significant  influence  on  the  computed  results  (Ref.  1) .  However, 
for  the  3-D  equations  with  large  axial  vorticity,  the  magnitude  of 
the  transverse  convective  velocities  can  be  of  the  same  order  as 
the  velocity  in  the  flow  direction.  In  this  situation,  the 
influence  of  the  model  for  the  diffusional  velocities  may  be  much 
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less  dominant  than  for  the  2-D  case.  Since  the  choice  of  model  for 
the  diffusional  velocities  has  a  first-order  impact  on  the 
computational  expense,  some  assessment  of  its  impact  on  the  rate  of 
mixing  in  flows  with  large  axial  vorticity  is  required. 

This  report  discusses  the  development  of  a  3-D  mixing  code  (TRIMIX) 
to  predict  the  influence  of  axial  vorticity  and  diffusional 
velocity  model  on  the  rate  of  mixing  between  the  primary  flow  and 
an  injected  secondary  flow.  The  parabolic  flow  approximation  will 
be  used  to  minimize  computation  time  and  the  initial  distribution 
of  axial  vorticity  will  be  specified.  Therefore,  the  calculation 
will  not  be  exact  since  the  initic.„  vorticity  distribution  will  not 
be  known  exactly.  However,  the  experimental  data  base  on  sidewall 
injection  is  sufficient  to  make  a  reasonable  estimate  of  the 
initial  vorticity.  The  resulting  code  is  expected  to  run  between 
one  and  two  orders  of  magnitude  faster  than  a  Navier-Stokes 
calculation  (3  to  5  h  of  CRAY  time  for  a  chemically  reacting  flow)  . 
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SECTION  2 
ANALYSIS 


A  typical  chemical  laser  nozzle  is  shown  in  Figure  2.1  with 
sidewall  injection  from  both  walls  in  the  subsonic  flow  region. 
Also  shown  is  a  commonly  used  injection  pattern  with  two  rows  of 
staggered  holes.  The  number  of  rows,  hole  sizes  and  spacing  are 
chosen  to  provide  a  uniform  distribution  of  the  secondary  flow 
throughout  the  primary  flow.  A  rectangular  coordinate  system  will 
be  used  with  x  in  the  primary  flow  direction,  y  perpendicular  to 
the  nozzle  axis  and  z  in  the  direction  along  which  the  injection 
holes  are  spaced.  The  corresponding  velocity  components  are  u  in 
the  x  direction,  v  in  the  y  direction  and  w  in  the  z  direction. 
The  number  of  holes  in  the  z  direction  generally  ranges  from  20  to 
50.  Since  it  takes  at  least  20  mesh  points  in  the  z  direction  to 
resolve  the  flow  field  from  one  hole  to  the  next,  and  about  40  in 
the  y  direction  from  y  =  o  to  y  =  H,  it  is  obvious  that  carrying 
out  the  computation  over  the  entire  array  of  holes  would  be 
prohibitively  expensive.  Fortunately,  in  most  instances  the  hole 
pattern  in  the  z  direction  is  symmetric  such  that  the  computations 
can  be  carried  out  for  a  unit  cell  of  width  A  (see  Fig.  2.1). 
Therefore,  the  computational  domain  in  the  z  direction  is  from  z  = 
o  to  z  =  A(x)  .  Since  the  walls  which  bound  the  flow  in  the  z 
direction  are  generally  not  parallel  to  one  another,  A  must  vary  in 
the  x  direction  to  conserve  the  cross-sectional  flow  area.  The 
injection  pattern  from  the  top  wall  (y  =  H)  is  usually  staggered 
from  the  hole  pattern  at  y  =  -H  such  that  the  flow  is  not  symmetric 
about  y  =  o.  However,  to  conserve  computation  time,  in  this  report 
y  =  o  will  be  assumed  to  be  a  symmetry  plane  and  the  computational 
domain  in  the  y  direction  is  from  y  =  o  to  y  =  H(x) . 

2 . 1  EQUATIONS . 

To  satisfy  the  continuity  equation,  two  stream  functions  are 
defined  as  follows: 
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Figure  2.1.  Nozzle  geometry  and  injection  hole  pattern 
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The  remaining  velocity  component  v  is  defined  by: 


pv  =  Hr  +  ^  ( 2 • 2 ) 

The  axial  vorticity  (i.e.,  the  x  vorticity  component)  is  defined 
by: 


, ,  _  3w  3v 
3y  az 


(2.3) 


where  the  vorticity  is  defined  to  be  positive  for  clockwise 
rotation. 


Substitution  of  the  v  and  w  expressions  of  Equations  2 . 1  and  2 . 2 
into  Equation  2.3  provides  the  equation  for  the  stream  function  Y. 


a2v  aff 
ay2  az2 


i  d£_)  df_ 
fi  dy)  3y 


( 1  &_)  av 

kp  az  J  az 


=  -  pw  + 


(2-4) 


The  momentum  equation  in  the  x  direction  is  used  to  solve  for  the 
u  velocity  component  and  is  given  by: 


pV-Vu  =  V-  (pVu)  -  (P  (2.5) 

IP  s 

dx 

This  equation  is  parabolized  by  dropping  the  x  derivative  in  the 
divergence  of  pVu  and  by  replacing  the  pressure  gradient  term  by 
the  derivative  of  the  mean  pressure.  That  is,  the  pressure 
gradient  (P)  is  some  average  value  over  the  flow  cross  section  ar.d 
is  chosen  to  satisfy  global  continuity  in  the  form  of  Equation  2.6. 
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(2.6) 


H  A 

^  =  pudzdy  =  constant 
o  o 

Equations  2.5  and  2.6  are  two  coupled  equations  which  are  solved 
for  u  and  <P.  After  solving  for  u,  the  stream  function  ip  of 
Equation  2.1  is  evaluated  by  integration: 


y 

<P  =  -  \  L(pu)dy  (2.7) 

J  o 

The  transport  equation  for  axial  vorticity  for  a  compressible, 
variable  viscosity  fluid  is  derived  in  Appendix  A  and  is  given  as 
(Equation  A. 6) : 

pV-Vw  =  V-(pVw)  -  cj  ^  +  Q  (2.8) 

where  the  subscript  x  has  been  dropped  from  wx  and  Equation  2.1  has 
been  used  to  introduce  the  stream  function  ip.  The  quantity  Q  is 
defined  in  Equation  A.  6  of  Appendix  A  and  contains  the  terms 
involving  the  gradients  in  density  and  viscosity. 

The  species  equations  are  given  by: 

PV-VK  i  =  -  V-fpKjV^  +  wjL  (2.9) 

where  Kf  is  the  species  mass  fraction,  w,-  is  the  net  rate  of 
production  of  species  i  due  to  chemical  reactions  and  V;  is  the 
species  diffusional  velocity  vector  and  is  defined  by: 


V. 

i 


v . 

l 


+  w.  e_ 

y  i  z 


(2 . 10) 


where  axial  diffusion  is  ignored  due  to  the  parabolic  flow 
approximation . 
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As  indicated  earlier,  one  of  the  objectives  of  this  effort  is  to 
assess  the  influence  of  the  diffusional  velocity  model  on  the 
computed  results.  In  Reference  2,  it  was  demonstrated  that  the 
"effective"  binary  diffusion  model  gave  results  which  were  in 
excellent  agreement  with  the  exact  multicomponent  diffusion  model. 
Since  the  "effective"  binary  diffusion  model  requires  much  less 
computational  effort  than  solving  the  Stefan-Maxwell  equations  for 
the  multicomponent  model,  the  effective  binary  diffusion  model  is 
used  to  define  the  diffusional  velocities.  The  diffusional 
velocities  are  given  as: 


K.Mv 


i 


KiMwi 


N 

Di7y(K.M)  +  Ki  I  Dj Vy (KjM)  +  S.VyP/P 
J=1 

N 

DiVz(KiM)  +  Ki  ^  DjVz(KjM)  +  SiVzP/P 
J=1 


Si 


Wi 


Mi/M) 


Ki 


N 

I 


DjKj(l 


Mj/M) 


J=1 


(2.11) 


where  Vy  denotes  the  gradient  operator  in  the  y  direction  and  Vz  in 
the  z  direction.  M  is  the  mixture  molecular  weight.  Reference  3 
may  be  consulted  for  the  derivation  of  Equation  2.11.  The 
diffusion  coefficient  is  given  by: 


D. 


l 


(1 


(x./pD.  .  ) 


(2.12) 


where  x(  is  the  mole  fraction  of  species  i  and  the  DM  are  the 
binary  diffusion  coefficients.  Note  that  the  summation  of  Equation 
2.12  excludes  j  =  i. 


For  comparative  purposes,  the  code  also  contains  the  binary 
diffusion  model  with  the  diffusional  velocities  given  by: 
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As  indicated  in  Section  2,  the  planes  y  =  o,  z  =  o  and  z  =  A  are 
symmetry  planes.  Therefore,  at  y  =  o,  the  y  gradients  in  axial 
velocity,  temperature  and  species  vanish.  At  z  =  o  and  z  =  A,  the 
z  gradients  in  axial  velocity,  temperature  and  species  vanish.  At 
y  =  H,  there  may  be  either  a  wall  or  a  plane  of  symmetry.  If  a 
symmetry  plane  is  present  at  y  =  H,  then  the  y  gradients  in 
velocity,  temperature  and  species  are  zero.  For  a  wall  at  y  =  H, 
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the  no-slip  condition  (u  =  o)  is  used  for  velocity  and  the  wall 
temperature  is  specified. 

Heterogeneous  wall  reactions  are  accounted  for  in  the  species 
boundary  conditions  at  y  =  H  when  a  wall  is  present.  The  boundary 
condition  is: 

y  =  H:  (for  wall) 

-  PwCos9  Ki«i  =  ji 

Tan0  =  d^  (2.16) 

The  net  rate  of  species  production  per  unit  time  and  surface  area 
is  denoted  by  jf  and  is  given  as: 


J 


i 


(  ) 

V2 

- 

(RT/277M.^ 

*lKi  +  *2_ 

(2.17) 


where  the  specific  expressions  for  (j>:  and  <pz  depend  on  the 
particular  wall  reactions  involving  the  ith  species.  For  atom  wall 
recombination  reactions,  the  surface  reaction  is: 


y 

A  +  ..  +  wall  - ►  A2  +  wall 


(2.18) 


For  reaction  Equation  2.18,  the  coefficients  are: 


i  =  A:  <p±  =  -  y 

i  =  a2:  <P1  =  -  yV 2 


*2  "  ^KA. 


Vka 


2j 


eq 


VKA. 


eq 


<P2  =  ka 


(2.19) 


where  y  is  the  recombination  coefficient.  Note  that  usually  the 
wall  temperature  is  low  enough  such  that  the  equilibrium  ratio  is 
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vanishing  small,  in  which  case  #2(i  =  A)  and  ^(i  =  A2)  may  be 
zeroed  with  negligible  error. 

For  surface  deactivation  of  an  excited  species  A*,  the  reaction  is 
represented  by: 

y 

A*  +  wall  - ►  A  +  wall 


i  =  A*: 

<P1  =  ~  Y 

<t>2  =  0.0 

i  =  A: 

(p±  =  0.0 

<t>2  =  yK(A*) 

(2.20) 

The 

boundary  conditions 

for  the  axial  vorticity 

and  stream 

functions  are  obtained  from  their  definitions  and  a  knowledge  of 

the 

behavior  of 

the  v  and 

w  velocity  components  on  the 

boundaries . 

The 

v  and  w  velocity  components  on  the  boundaries  are 

given  by: 

y  =  o:  v 

=  o 

3w 

dy  ~  ° 

y  =  H: 

v-  =  u 

dx 

w  =  o  (wall) 

9w 

=  o  (symmetry  plane) 

3v 

z  ~  o:  dz 

=  o 

w  =  o 

z  =  4:  §f 

=  o 

dA 

W  =  U  3 — 

dx 

(2.21) 

At  the  symmetry  planes  y  =  o  and  z  =  o,  the  axial  vorticity  is 
zero.  At  z  =  A,  differentiating  the  boundary  condition  on  w  with 
respect  to  y  gives: 


z  =  A: 


,,  _  dA  <9u 
dx  <9y 


(2.22) 


At  y  =  H,  the  axial  vorticity  is  zero  if  a  symmetry  plane  is 
present.  If  a  wall  is  present  at  y  =  H,  then  v  =  w  =  o  and  the 
vorticity  is: 
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y  =  H: 


w 


-i  a2¥ 


(2.23) 


p  3y2 

The  boundary  condition  on  ip  at  y  =  o  must  be  consistent  with  the 
requirement  that  v  =  o  at  the  symmetry  plane.  Therefore,  ip  =  o  at 
y  =  o.  This  condition  is  incorporated  into  the  expression  of 
Equation  2.7.  The  derivatives  of  ¥  with  respect  to  y  and  z  are 
evaluated  on  the  boundaries  using  Equations  2.1,  2.2  and  2.21.  By 
integrating  these  derivatives,  the  ¥  values  on  the  boundaries  can 
be  determined  to  within  an  arbitrary  constant.  Since  only  the 
derivatives  of  ¥  are  required  to  define  v  and  w,  the  value  assigned 
to  the  constant  is  irrelevant.  Therefore,  we  set  ¥  =  o  at  y  =  o, 
z  =  o  to  obtain  the  following: 


y  =  o:  ¥  =  o 


y  - 


H: 


z  =  o: 


z  =  A: 


¥ 


dA 

dx 


y 

pudy 

o 


2.3  INITIAL  CONDITIONS. 


(2.24) 


Figure  2.2  shows  the  injected  jets  corresponding  to  the  hole 
injection  patterns  of  Figure  2.1.  In  the  immediate  region  of  the 
injection  process  (between  stations  0  and  D)  the  problem  is 
elliptic  and  requires  a  Navier-Stokes  solution  to  define  the  flow 
field.  The  approach  in  the  TRIMIX  code  is  to  start  the  calculation 
at  some  small  distance  downstream  (station  D)  of  the  jet  injection 
location  where  the  primary  flow  has  bent  the  jets  approximately 
parallel  to  the  flow.  Since  the  jets  are  still  essentially  unmixed 
at  station  D,  this  seems  to  be  a  reasonable  way  to  carry  out  the  3- 
D  mixing  calculation  while  avoiding  the  complexity  of  the  exact 
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Navier-Stokes  calculation.  The  difficulty  with  this  approach  is 
that  it  requires  the  specification  of  the  initial  conditions  at 
station  D.  In  Appendix  B,  the  initial  conditions  for  the 
thermodynamic  state  variables,  the  axial  velocity  (u)  and  the  area 
occupied  by  each  jet  are  obtained  by  satisfying  the  global 
conservation  equations  with  some  closure  assumptions.  To  obtain 
the  initial  axial  vorticity  and  transverse  velocity  components  (v 
and  w) ,  proceed  as  follows. 

Following  Fearn  and  Weston  (Ref.  4) ,  we  adopt  the  diffuse  vortex 
model  where  the  distribution  of  vorticity  is  assumed  to  be  Gaussian 
within  each  vortex.  It  is  also  assumed  that  the  jets  are  spaced 
far  enough  apart  such  that  the  vorticity  field  of  each  jet  is 
independent  from  that  of  its  neighbors.  Thus,  the  vorticity 
distribution  is  taken  to  be  a  linear  combination  of  the  vorticity 
field  of  each  jet.  This  gives: 


w 


l 


i=l 


uj  ^e 


(  )2  (  ) 

=  y-YiJ  +  l*-«ij 


(2.25) 


where  r  is  the  radial  distance  from  the  center  of  each  vortex  with 
coordinates  y(  and  z-.  For  the  injection  pattern  shown  in  Figure 
2.1,  there  are  three  vortices  located  in  the  unit  computational 
cell.  The  vortex  locations  are  shown  in  Figure  2.3. 


The  strength  of  each  vortex  (f)  is  defined  by: 


r . 

l 


2  n 


rdr  =  nui/pi 


o 


(2.26) 


The  diffusion  constant  /?  is  defined  in  terms  of  the  vortex  core 
size  by  Equation  2.6  of  Reference  4  and  is  reproduced  here  as: 
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(1^121) 


P  = 


where  the  P  of  this  report  is  p  of  Fearn  and  Weston. 


(2.27) 


The  vortex  core  size  is  obtained  from  Figure  9  of  Reference  4  and 
is  curve  fit  as: 


0.11  |  +  0.60 


(2.28) 


where  S  is  distance  measured  along  the  vortex  curve  from  the 
injection  location  to  x  =  o.  The  distance  S  is  defined  in  Equation 
2.33. 


The  vorticity  at  the  center,  is  obtained  from  Equation  2.26  for 
a  given  value  of  the  vortex  strength  I",- .  The  data  of  Fearn  and 
Weston  (Table  2  of  Ref.  4)  for  the  initial  vortex  strength  near  the 
jet  orifice  is  plotted  in  Figure  2.4  for  normal  injection  into  a 
subsonic  flow.  Also  shown  is  the  maximum  computed  vorticity  from 
the  analytical  model  of  Karagozian  (Ref.  5).  The  data  are 
represented  by: 


r  =  1.4  DV  R 
00 


R  = 


(2.29) 


where  D  is  the  jet  orifice  diameter,  subscript  J  denotes  jet 
orifice  values  and  subscript  00  denotes  free  stream  values  (i.e., 
station  O  in  Fig.  2.2).  Note  that  and  I",-  are  defined  to  be 
positive  for  clockwise  rotation  and  negative  for  counterclockwise 
rotation. 
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10.0 


Figure  2.4. 


Initial  vortex  strength  for  normal  injection  into 
subsonic  flow. 
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The  location  of  each  vortex  center  is  estimated  from  the  jet 
centerline  and  vortex  center  trajectory  data  of  Fearn  and  Weston 
(Ref.  4) .  The  general  form  of  the  trajectory  equation  is  given  by 
the  following  expression  where  the  a,  b,  and  c  coefficients  are 
different  for  the  jet  centerline  and  vortex  center. 

h/D  =  aRb  (x/D)c  (2.30) 

where  h  is  distance  from  the  wall  to  either  the  jet  centerline  or 
the  vortex  center.  The  coefficients  are  given  for  the  jet 
centerline  as:  a  =  0.9772,  b  =  0.9113,  c  =  0.3346  and  for  the 
vortex  center  as:  a  =  0.3473,  b  =  1.127  and  c  =  0.4291. 


The  angle  of  inclination  (0)  of  the  jet  centerline  to  the  wall  is 
obtained  from  the  slope  of  Equation  2.30.  At  x  =  o,  the  largest 
angle  is  d2  (for  jet  2,  see  Fig.  2.2).  Thus,  for  a  specified  value 
of  02,  the  corresponding  x2  distance  is  given  by: 


x2  "  d2 (° • 327R0’ 911/Tan02 j 1  * 5 


(2.31) 


The  x1  location  is  given  by: 


Xx  =  x2  +  L  (2.32) 

For  a  given  value  of  9Z  (typically  20  deg),  Equations  2.31  and  2.32 
define  x.,/0.,  and  x2/D2  and  Equation  2.3  0  gives  the  values  of  and 
h2  for  the  vortex  center  and  the  jet  centerline  location  at  x  =  o. 


The  distance  along  the  vortex  center  curve  is  given  by: 


S  = 


Cl  ,  fdhl 

2- 

L1  U*J 

1/2 


dx 


fj  =  O.USR1-127/!^))0-571 


(2.33) 


where  x,  defines  S1  and  x2  defines  S2. 
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The  vortex  half-spacing  (d,,  d2,  d3)  of  Figure  2.3  is  obtained  from 
Figure  7  of  Reference  4  and  is  given  in  Figure  2.5.  The  x 
locations  for  the  data  of  Figure  2.5  were  determined  for  02  =  20 
deg. 

The  actual  jet  shape  is  that  of  a  horseshoe  with  the  jet  centerline 
at  the  closed  end  and  the  vortices  located  in  the  legs  at  the  open 
end.  Although  there  is  no  intrinsic  difficulty  in  using  a  fine 
mesh  to  resolve  the  horseshoe  geometry,  the  coding  to  define  the 
mesh  becomes  quite  complicated  and  the  computational  cost  is 
increased  significantly.  Therefore,  the  jet  geometry  has  been 
approximated  by  an  elliptical  cross  section  with  the  vortex  located 
on  the  major  axis  of  the  ellipse.  The  actual  jet  geometry  used  in 
the  calculation  is  shown  in  Figure  2.6  where  the  elliptical  cross 
section  has  been  truncated  into  a  rectangle  to  simplify  the  mesh 
generation.  The  unit  computational  cell  in  Figure  2.6  corresponds 
to  that  shown  in  Figures  2.1  and  2.2  and  contains  one-half  of  a 
large  jet  (JET1)  and  two  halves  of  the  small  jet  ( JET2) .  The  jet 
areas  are  defined  in  Appendix  B  and  the  location  is  defined  by 
Equations  2.30  through  2.32.  To  close  the  geometry,  the  ratio  of 
jet  major  to  minor  axis  (€)  must  be  specified.  Experimental  data 
suggest  values  of  6  in  the  range  of  one  to  five.  However,  there 
are  additional  constraints  on  G  dictated  by  the  hole  spacing  and 
the  vortex  locations.  That  is,  the  extent  of  the  jets  in  the  z 
direction  must  be  large  enough  to  incorporate  the  vortices  inside 
the  jet  and  small  enough  such  that  the  two  half- jets  do  not  touch. 
With  these  constraints,  a  value  of  G  =  2.0  seems  to  be  a  reasonable 
choice  and  this  value  has  been  used  for  all  the  calculations  of 
this  report. 

The  experimental  information  given  above  should  provide  a 
reasonable  estimate  for  the  initial  jet  vorticity.  However,  since 
the  experimental  data  are  for  a  single  pressure  matched  jet 
exhausting  into  an  unbounded  flow,  it  is  important  to  emphasize 
that  this  is  only  an  estimate.  In  particular,  the  symmetry 
conditions  imposed  by  the  injection  hole  spacing  and  locations  are 
not  present  in  the  experimental  data  base. 
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Figure  2.5.  Vortex  half  spacing. 
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JET  1 


SECTION  3 

NUMERICAL  ANALYSIS 


The  equations  of  Section  2  are  defined  in  terms  of  the  physical 
coordinates  x,  y,  and  z.  The  actual  calculations  are  carried  out 
in  the  transformed  r]  and  £  coordinate  system.  The  transformed 
coordinates  are  defined  by: 

V  =  Y/H(x)  £  =  z/A(x)  (3.1) 

where  a  nonuniform  cartesian  mesh  is  used  in  the  rj -£  coordinate 
system. 

The  equations  are  differenced  using  a  first-order  accurate  two- 
point  backward  difference  for  the  x  derivatives  and  second-order 
accurate  three-point  central  difference  for  the  77  and 
differences. 

After  differencing,  all  the  linearized  equations  have  the  following 
form: 


AVi,k  +  BDj,k  +  cuj-i,k  +  DUj,k+i +  EUj,k-i +  ^  -  F  <3-2> 


where  the  coefficients  A,  B,  C,  D,  E  and  F  are  all  known  functions 
of  j  and  k.  The  index  j  is  in  the  r]  direction  and  the  k  index  is 
in  the  £  direction.  The  term  with  <P  denotes  the  pressure  gradient 
in  the  axial  momentum  equation  (Eq.  2.5)  and  is  not  present  for  the 
other  equations.  Therefore,  w  =  -1.0  for  the  momentum  equation  and 
w  =  o  for  all  the  other  equations.  The  variable  U  is  solved  for 
using  an  implicit  formulation  where  the  block  tridiagonal  approach 
is  used  to  solve  the  linear  system  of  equations.  An  iterative 
approach  is  used  where  the  coefficient  arrays  are  updated  after 
each  iteration  and  Equation  3.2  is  solved  repeatedly  until  some 
error  criterion  is  satisfied.  Then  a  new  step  is  taken  in  the  x 
direction  and  the  process  is  repeated. 
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When  the  momentum  equation  is  solved,  the  pressure  gradient  (P  is 
solved  simultaneously  with  the  axial  velocity  by  satisfying 
Equation  2.6  in  the  form: 


1  1 

0=|  |  0U  d^d£ 

o  o 

<p  =  m/  (HA)  (3.3) 

where  9  is  a  given  function  and  U  denotes  the  axial  velocity.  For 
incompressible  or  low-speed  flow  where  the  density  varies  slowly 
with  pressure,  9  is  the  density  p  evaluated  at  the  previous  iterate 
level.  For  supersonic  flow,  where  the  density  varies  strongly  with 
pressure.  Equation  2.6  is  written  as: 


0U  ctyd£ 


9  =  M/RT  (3.4) 

where  T  is  the  temperature  and  M  is  the  mixture  molecular  weight. 
The  pressure  P  is  related  to  the  pressure  gradient  (P  through  the 
use  of  a  backward  difference: 


<P  =  (P  -  PL)/dx  (3.5) 
where  PL  is  the  pressure  at  the  previous  x  station. 

Thus,  for  incompressible  or  subsonic  flow,  Equations  3.2  and  3.3 
are  solved  simultaneously  for  U  and  <?.  For  supersonic  flow. 
Equations  3.2,  3.4  and  3.5  are  solved  simultaneously  for  U,  (P  and 
P. 


The  solution  procedure  is  to  use  some  quadrature  to  express  the 
integral  of  Equation  3.3  or  3.4  as  a  linear  sum  which  can  be  solved 
within  the  framework  of  the  block  tridiagonal  approach.  This 
approach  was  used  by  Crowell  (Ref.  1)  and  Truman  (Ref.  6)  for  2-D 
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flows.  The  extension  to  3-D  flows  is  straightforward  but  tedious 
and  the  details  are  omitted  here.  The  difference  between  the 
solution  using  the  subsonic  or  supersonic  formulations  is  that  the 
supersonic  formulation  results  in  a  quadratic  for  the  pressure, 
whereas  the  subsonic  or  constant  density  formulation  yields  only  a 
single  solution  branch.  The  two  solution  branches  obtained  from 
the  use  of  Equation  3.4  are  both  physically  valid  and  represent 
subsonic  and  supersonic  solutions.  For  supersonic  flow,  the  lower 
pressure  is  always  used.  For  low  subsonic  flows,  the  higher 
pressure  is  always  used.  However,  for  subsonic  flow  near  Mach  1, 
the  correct  pressure  can  be  on  either  solution  branch  and  some  care 
is  required  in  choosing  the  correct  value. 

The  solution  procedure  for  Equation  3.2  using  the  block  tridiagonal 
approach  marches  in  the  j  direction  from  J  =  1  to  J  =  JMAX  and 
solves  a  full  matrix  for  the  KMAX  values  of  the  recursion 
coefficients  at  each  J  value.  The  solution  for  U  is  obtained  by 
oack  substitution  from  J  =  JMAX  to  J  =  1.  The  computational  effort 
is  proportional  to  JMAX (KMAX)3.  Therefore,  it  is  extremely 
important  to  choose  the  J  direction  as  the  direction  with  the 
largest  number  of  mesh  points.  For  the  case  shown  in  Figures  2.1 
through  2.3,  the  H  dimension  is  nearly  four  times  the  A  dimension 
and  there  are  about  three  times  as  many  mesh  points  in  the  y 
direction  as  in  the  z  direction.  Therefore,  it  is  critical  that 
the  J  dimension  be  in  the  y  direction.  Choosing  the  J  dimension  in 
the  z  direction  would  increase  computational  cost  by  a  factor  of  9. 
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SECTION  4 

RESULTS  AND  DISCUSSION 


The  case  of  most  interest  to  WL/ARDI  is  the  modeling  of  the  flow  in 
the  RotoCOIL  device  where  the  iodine  injection  takes  place  in  a 
subsonic  channel  and  the  flow  is  expanded  through  a  throat  to 
achieve  supersonic  cavity  flow.  In  principle,  the  TRIMIX  code  can 
model  this  flow  geometry  by  iterating  upon  the  upstream  pressure 
until  sonic  flow  in  the  throat  is  achieved.  However,  this 
constitutes  a  classical  problem  in  fluid  dynamics  where  a  saddle 
point  singularity  exists  at  Mach  1  and  can  only  be  integrated  with 
great  difficulty,  using  a  marching  technique.  The  practical 
difficulty  is  that  the  solution  must  be  extremely  accurate  to 
integrate  smoothly  through  the  throat.  This  translates  into  a 
requirement  that  the  initial  pressure  must  be  iterated  until  it  is 
known  to  within  a  tolerance  of  0.01  torr.  Typically,  this  would 
take  10  or  12  iterations,  each  one  of  which  would  consume  3  to  4  h 
of  computer  (CRAY  1)  time.  Although  the  solution  obtained  in  this 
manner  is  exact,  it  is  computationally  expensive.  One  way  around 
this  problem  is  to  specify  the  pressure  distribution  and  accept  the 
associated  error  in  the  global  mass  flux.  That  is,  if  (P  is 
specified  in  the  momentum  equation,  then  there  is  no  way  to  ensure 
that  Equation  2.6  is  satisfied.  That  is,  an  error  in  pressure  will 
result  in  an  error  in  the  global  mass  flux.  Therefore,  to  limit 
the  mass  flux  error  to  some  acceptable  level,  a  reasonably  accurate 
pressure  distribution  must  be  supplied  to  the  calculation.  Since 
the  pressure  is  determined  primarily  by  the  area  schedule  and  the 
loss  in  momentum  due  to  wall  shear,  it  seems  plausible  to  assume 
that  a  one-dimensional  (1-D)  premixed  calculation  could  be  used  to 
obtain  the  pressure  distribution.  The  1-D  premixed  calculation 
must  also  iterate  on  the  initial  pressure  to  achieve  sonic  flow  at 
the  throat  but  the  computational  cost  of  this  calculation  is  less 
than  a  half-hour  of  CRAY  1  time.  This  approach  has  been  used  with 
the  2-D  JETMIX  code  for  mixing  calculations  and  has  worked  quite 
well,  with  the  error  in  mass  flux  oscillating  between  5  and  10 
percent  over  all  the  streamwise  locations.  Since  the  mass  flux 
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error  is  distributed  over  all  the  species,  it  is  believed  that  the 
error  in  the  species  mass  fractions  is  less  than  the  error  in  the 
global  mass  flux.  This  approach  should  also  work  for  the  3-D 
mixing  problem  and  therefore  the  TRIMIX  code  contains  an  option  to 
either  calculate  the  pressure  by  satisfying  global  continuity  or 
specify  the  pressure. 

Since  it  is  believed  that  most  of  the  mixing  and  I2  dissociation 
takes  place  in  the  subsonic  flow  upstream  of  the  throat,  the 
calculations  will  concentrate  on  this  region.  To  avoid  the 
pressure  iteration  described  above,  the  throat  is  removed  and  the 
calculations  are  done  for  a  constant  area  channel  where  the  initial 
conditions  are  those  of  the  RotoCOIL  nominal  operating  point.  The 
conditions  upstream  of  the  iodine  injection  and  the  geometry  are 
given  as: 


p 

s 

34.71  torr 

T 

= 

275  K 

Y 

= 

0.37 

U 

= 

0.93 

P(H20) 

= 

0.897  torr 

x(He) 

SB 

4.54  moles/s 

x(02) 

= 

1.385 

x(Cl2) 

- 

0.107 

x(H20) 

= 

0.16 

m 

= 

0.023186  g/s 

H 

= 

0.37338  cm 

A 

SB 

0.09726  cm 

Note  that  the  molar  flow  rates  are  those  for  the  entire  nozzle  bank 
whereas  the  mass  flux  (m)  is  for  the  unit  computational  cell 
defined  by  H  and  A.  The  I2/02  molar  flow  ratio  is  0.01487  and  the 
injected  molar  ratio  of  He/I2  is  43.835.  It  is  assumed  that  the 
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injected  flow  is  choked  at  the  orifice  hole  to  give  a  value  of  the 
square  root  of  the  momentum  ratio  of  R  -  2.763.  The  injected  I2 
and  He  has  a  total  temperature  of  366  K.  The  corresponding 
information  from  Section  2.3  is  given  as: 


t/dv. 

= 

3.87 

a 

0.3865 

Vz 

= 

0.747 

*pi 

a 

0.2612 

<2 

a 

0.5223 

where  the  r]  coordinates  are  the  jet  centerline  locations  for  the 
indicated  jets  and  the  coordinates  are  the  vortex  center 
locations  (see  Figs.  2.3  and  2.6).  The  initial  jet  and  vortex 
locations  are  shown  in  the  transformed  plane  in  Figure  2.6.  For 
all  the  calculations  of  this  report,  10  species  (I2,  He,  I,  I2*,  I*, 
H20,  O^A),  02 ( )  ,  02(3I)  ,  Cl2)  and  the  21  reactions  of  the  reduced 
WL  oxygen/iodine  kinetics  package  (Table  2.1-VI  of  Ref.  7)  were 
used. 

The  calculations  used  48  mesh  points  in  the  r?  direction  and  16  mesh 
points  in  the  £  direction.  The  step  size  in  the  x  direction 
started  at  0.01  cm  and  was  ramped  up  to  0.03  cm.  The  calculation 
was  run  out  to  1.8  cm  and  a  total  of  66  steps  were  taken  in  the  x 
direction.  The  computation  time  was  3.89  h  on  the  CRAY  1  machine. 

The  initial  distribution  for  the  axial  vorticity  and  the  velocity 
components  are  contained  in  Figure  4.1.  Note  that  the  vorticity  is 
large  enough  to  induce  secondary  flow  velocities  of  the  same 
magnitude  as  the  axial  velocity.  The  convection  and  diffusion  of 
the  iodine  jet  can  be  seen  in  Figure  4.2  with  increasing  distance 
from  the  initial  x  location.  The  buildup  of  the  gain  is  shown  in 
Figure  4.3.  The  profiles  at  the  x  =  1.8  cm  station  are  contained 
in  Figure  4.4. 
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VORTICITY  MAP 


Figure  4.4.  Distributions  at  x  =  1.809  cm. 


TEMPERATURE  MAP 


4.1  COMPARISON  OF  3-D  MIXING  AND  1-D  PREMIXED  CALCULATIONS. 


The  premixed  assumption  is  frequently  used  to  carry  out  performance 
calculations  for  COIL  laser  devices  (see  Ref.  8  for  example) .  In 
this  manner,  the  mixing  calculation  is  eliminated  by  assuming 
instantaneous  mixing  between  the  primary  and  secondary  flows.  This 
reduces  the  3-D  chemically  reacting  flow  problem  to  a  1-D  problem 
with  a  corresponding  reduction  in  computer  cost  by  a  factor  of  100. 
The  TRIMIX  code  contains  an  option  to  premix  the  reactants  at  the 
initial  station  and  carry  out  the  premixed  calculation  to  provide 
a  comparison  between  the  premixed  assumption  and  the  exact  3-D 
results.  To  compare  the  results,  average  values  are  defined  for 
the  gain  and  species  as  follows: 
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(4.3) 


where  Kf  and  or  are  the  species  mass  fractions  and  gain  and  F  is  the 
fraction  of  molecular  iodine  which  has  been  dissociated.  The  bar 
superscript  denotes  values  which  are  averaged  over  the  entire 
computational  domain  at  each  x  station.  The  average  gain  in  the 
direction  of  the  optical  axis  (y-direction)  determines  the 
radiative  flux  and  is  defined  by: 
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where  the  gain  calculation  is  described  in  Appendix  C. 


(4.4) 
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The  comparison  is  contained  in  Figures  4.5  and  4.6.  Figure  4.5 
shows  the  fraction  of  I2  dissociated  (F  of  Eq.  4.3)  and  the  average 
gain  (tf)  as  a  function  of  flow  distance  and  Figure  2 . 12  shows  the 
optical  axis  gain  (Eq.  4.3).  The  premixed  assumption  drastically 
underpredicts  both  the  amount  of  iodine  dissociation  and  the  small 
signal  gain.  To  the  author's  knowledge,  this  is  the  first 
comparison  of  the  premixed  assumption  with  an  exact  3-D  mixing 
calculation  for  an  oxygen/iodine  laser. 

The  opinion  is  often  expressed  that  while  the  premixed  calculation 

is  not  exact,  it  does  provide  the  upper  bound  (i.e.,  most 

optimistic)  laser  performance.  The  results  of  Figures  4.5  and  4.6 

clearly  dispute  this  assertion.  In  this  case,  the  premixed 

assumption  substantially  underpredicts  the  performance.  This 

difficulty  with  the  premixed  assumption  has  been  realized  for  some 

time  (see  the  discussion  on  pages  25-28  of  Ref.  8) .  The  approach 

used  in  Reference  8  to  overcome  this  limitation  was  to  artificially 

★ 

increase  the  rate  constant  for  the  I2  formation  reaction.  That  is, 
the  first  step  in  the  I2  dissociation  chain  is  the  following 
reaction. 

I2  +  02(14)  - ►  I2*  +  02(3I) 

K  =  7.0x10  15  cm3/molecule-s  (4.5) 

where  the  indicated  rate  is  the  "standard"  value  from  Reference  7. 
Both  the  3-D  mixing  calculation  and  the  premixed  calculations  of 
Figures  4.5  and  4.6  used  the  rate  constant  of  Equation  4.5. 

If  the  I2  dissociation  fraction  is  known  at  some  x  location,  then 
the  rate  constant  of  Equation  4.5  can  be  chosen  to  reproduce  the 
known  dissociation.  From  the  3-D  calculation,  the  dissociation 
fraction  is  seen  to  be  about  44  percent  at  x  =  1.83  cm.  By 
increasing  the  rate  constant  of  Equation  4.5  by  a  factor  of  5,  the 
premixed  calculation  will  yield  the  same  result  at  x  =  1.83  cm. 
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Figure  4.5. 


Comparison  of  3-D  mixing  and  1-D  premixed 
calculations  for  I2  dissociation  and  average  gain 
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Figure  4.6. 


Comparison  of  3-D  mixing  and  l-D  premixed 
calculations  for  optical  axis  gain. 


This  is  shown  in  Figure  4.7.  The  corresponding  small  signal  gain 
and  yield  are  compared  with  the  average  values  from  the  3-D 
calculation  in  Figure  4.8,  and  the  agreement  is  quite  good.  In 
other  words,  if  the  rate  constant  of  Equation  4.5  is  increased  to 
reproduce  the  known  level  of  I2  dissociation,  then  the  premixed 
calculation  will  give  the  correct  average  gain  and  yield.  The 
comparison  of  Figures  4.7  and  4.8  seems  to  provide  the 
justification  for  the  use  of  the  premixed  approximation  with  the 
approach  of  Reference  8. 

4.2  COMPARISON  OF  DIFFUSION  MODELS. 

The  binary  diffusion  model  (Eq.  2.13)  is  frequently  used  in  mixing 
calculations  because  of  its  simplicity.  One  of  the  objectives  of 
this  report  was  to  make  a  comparison  between  the  binary  diffusion 
and  "effective”  binary  diffusion  (Eq.  2.11)  models.  This 
comparison  is  contained  in  Figures  4.9  and  4.10.  In  Figure  4.9, 
the  iodine  dissociation  and  optical  axis  gain  are  given  for  the  two 
models  where  a  Schmidt  number  of  unity  has  been  used  in  the  binary 
diffusion  model.  These  are  significant  differences  between  the  two 
models,  particularly  with  respect  to  the  iodine  dissociation 
fraction.  The  gains  profiles  at  x  =  1.8  cm  are  contained  in  Figure 
4.10,  which  shows  that  the  profile  for  the  binary  diffusion  model 
is  much  smoother  than  the  result  for  the  "effective"  binary 
diffusion  model. 

Although  the  differences  between  the  two  models  are  not  as  great  as 
those  reported  in  Reference  1,  they  are  still  sufficiently  large  to 
preclude  the  use  of  the  binary  diffusion  model.  There  undoubtedly 
exists  some  value  of  the  Schmidt  number  which  would  result  in 
closer  agreement  between  the  two  models.  The  difficulty  is  that 
the  "correct"  value  is  not  known  a  priori.  The  use  of  the 
effective  binary  diffusion  model  imposed  only  a  5  percent 
computational  penalty  for  the  calculations  of  Figures  4.9  and  4.10. 
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Figure  4.7 


Influence  of  I2  formation  rate  on  iodine  dissociation 
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Figure  4.8.  Influence  of  I2  formation  rate  on  gain  and  yield 
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Figure  4.9.  Influence  of  diffusion  model  on  gain  and  iodine 
dissociation. 
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Figure  4.10.  Influence  of  diffusion  model  on  gain  (x  =  1.80  cm). 
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Since  the  computational  penalty  is  small  and  there  is  no  ambiguity 
associated  with  choosing  a  value  of  the  Schmidt  number,  it  is 
recommended  that  the  "effective1*  binary  diffusion  model  be  used  in 
3-D  mixing  calculations. 

4.3  SENSITIVITY  OF  RESULTS  TO  UNCERTAINTIES  IN  INITIAL  CONDITIONS. 

As  indicated  in  Section  2.3,  the  experimental  data  base  upon  which 
the  initial  conditions  are  based  is  for  a  single  jet  issuing  into 
an  unbounded  flow.  Since  the  flow  geometry  for  the  present  problem 
consists  of  multiple  jets  being  injected  into  a  finite  flow,  the 
initial  conditions  of  Section  2 . 3  should  be  viewed  as  a  reasonable 
approximation  rather  than  as  an  exact  description  of  the  flow  at  x 
=  o.  In  order  to  assess  the  impact  of  uncertainties  in  the  initial 
conditions  upon  the  computed  results,  the  parameters  of  Section  2.3 
are  varied  over  a  plausible  range  and  the  calculations  are  compared 
with  the  nominal  case. 

The  (i  coefficients  of  Equation  2.25  determine  the  peak  value  and 
the  spread  of  the  vorticity.  As  P  increases,  for  a  fixed  value  of 
the  vortex  strength  V,  the  peak  vorticity  increases  and  the 
vorticity  is  more  concentrated  about  the  vortex  center.  The  p 
coefficients  are  defined  in  terms  of  the  vortex  core  size  in 
Equations  2.27  and  2.28.  If  the  vortex  core  radius  was  actually 
half  the  value  given  by  Equation  2.28,  the  P  coefficients  would  be 
four  times  as  large  as  the  nominal  values  and  the  peak  vorticity 
would  also  be  increased  by  a  factor  of  4 .  The  computed  results  for 
this  case  are  denoted  by  the  circle  symbols  in  Figure  4.11.  It  is 
seen  that  this  large  variation  in  p  has  only  a  minor  impact  on  the 
gain  and  I2  dissociation. 

There  is  also  some  uncertainty  in  the  vortex  half-spacing  (the  d 
coefficients  of  Fig.  2.3)  given  by  Figure  2.5.  The  nominal  values 
(given  by  £p1  and  ^  of  Eq.  4.2)  are  reduced  to  move  the  vortex 
locations  closer  to  the  center  of  the  jet.  The  nominal  vortex 
centers  are  shown  in  Figure  2.6.  Moving  the  vortex  center  to  the 
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midpoints  of  the  jets  gives  £p1  =  0.15625  and  ^  =  0.3125.  The 
calculated  results  for  this  variation  are  denoted  by  the  triangle 
symbols  in  Figure  4.11.  The  influence  on  the  results  is  seen  to  be 
negligible. 

The  variation  on  the  vortex  strength  is  contained  in  Figure  4.12 
where  the  nominal  value  (T/DV.  =  3.87)  has  been  approximately 
doubled  and  halved.  The  influence  on  the  gain  is  fairly  small  as 
is  the  influence  on  I2  dissociation  for  the  lower  value.  The 
effect  on  I2  dissociation  is  larger  for  T/DV.  =  8.0  (about  33 
percent  at  x  =  1.8  cm).  Even  so,  this  is  a  relatively  small 
variation  considering  the  magnitude  of  the  variation  in  the  vortex 
strength. 

Of  all  the  experimental  jet  data,  the  jet  trajectory  information  is 
the  most  accurate.  If  the  predicted  jet  penetration  is  close  to 
the  nozzle  centerline  ( r)  =  o)  ,  the  symmetry  condition  at  that 
location  would  probably  influence  the  jet  trajectory.  However,  for 
the  conditions  used  here  (Equation  4.1)  the  jet  closest  to  the 
centerline  (JET1)  is  still  well  away  from  it.  Therefore,  it  is 
believed  that  the  values  used  for  the  initial  jet  locations 
(Equation  2.30)  are  probably  accurate.  However,  to  see  what  impact 
the  symmetry  boundary  condition  at  T)  =  o  would  have,  a  case  was  run 
where  the  jet  centerlines  were  moved  towards  r]  =  o  such  that  the 
large  jet  (JET1)  was  just  touching  the  7=0  boundary  (r/1  =  0.0814 
and  T]z  =  0.44153).  The  results  indicated  negligible  influence  on 
the  I2  dissociation  but  the  gain  was  strongly  affected  with  the 
average  gain  being  reduced  from  0.0043/cm  to  0.0026/cm  at  x  =  1.8 
cm.  It  is  apparent,  therefore,  that  for  extreme  variations  in  jet 
location,  the  calculated  results  will  be  strongly  influenced. 

It  is  believed  that  the  results  of  this  section  demonstrate  that 
the  computed  results  are  only  moderately  influenced  by  relatively 
large  variations  in  the  initial  vorticity  field.  Large  variations 
in  jet  location  can  influence  the  results,  but  there  is  little 
uncertainty  in  the  jet  location. 
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Figure  4.11.  Influence  of  initial  conditions  on  gain  and  iodine 
dissociation. 


4.4  ZERO  VORTICITY  RESULTS. 


Figure  4 . 12  showed  that  reducing  the  nominal  vorticity  by  a  factor 
of  2  had  only  a  small  influence  on  the  computed  results.  The 
obvious  question  is,  What  would  be  the  impact  of  zeroing  the 
vorticity?  If  \_ne  results  are  only  slightly  affected  when  the 
initial  vorticity  is  zeroed,  it  would  suggest  that  the  transverse 
velocity  components  could  be  ignored  and  the  problem  treated  as  a 
3-D  diffusion  problem  with  a  single  convective  velocity  component 
(the  axial  velocity) .  The  comparison  with  the  initial  vorticity 
zeroed  is  contained  in  Figure  4.13,  which  shows  that  the  gain  is 
strongly  affected  by  zeroing  the  vorticity.  It  is  concluded  that 
the  flow  problem  cannot  be  simplified  by  using  a  zero  vorticity 
approximation . 
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Figure  4.12. 


Influence  of  initial  vortex  strength  on  iodine 
dissociation  and  gain. 
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SECTION  5 

SUMMARY  AND  CONCLUSIONS 


A  3-D  fluid  dynamic  computer  code  (TRIMIX)  has  been  written  to 
analyze  the  mixing  between  a  primary  flow  and  a  secondary  flow 
injected  as  a  series  of  discrete  jets  from  the  walls  which  bound 
the  flow.  The  flow  problem  has  been  parabolized  by  starting  the 
calculation  just  downstream  of  the  jet  injection  location  where  the 
jets  have  been  turned  parallel  to  the  primary  flow.  The  initial 
conditions  at  this  location  have  been  defined  by  assuming  that  the 
primary  and  secondary  flows  are  pressure  matched  and  unmixed.  The 
equations  are  solved  using  a  vorticity-stream  function  formulation 
and  the  initial  distribution  of  axial  vorticity  is  obtained  from 
the  experimental  data  base  on  jet  injection. 

A  number  of  calculations  have  been  carried  out  in  a  constant  area 
flow  channel  which  has  been  chosen  to  represent  the  subsonic 
portion  of  a  RotoCOIL  laser  nozzle.  The  oxygen/iodine  set  of 
species  and  reactions  consisting  of  10  species  and  21  reactions  has 
been  used  for  the  calculation.  The  equations  were  solved  over  a 
rectangular  grid  containing  768  mesh  points,  and  the  required 
computation  time  on  the  CRAY  1  computer  was  about  2  h/cm  of 
distance  in  the  flow  direction. 

The  calculations  indicate  that  moderate  uncertainties  in  the 
initial  conditions  have  only  a  minor  influence  on  the  computed 
results.  This  finding  provides  the  justification  for  the  approach 
used  in  the  TRIMIX  code;  namely,  that  the  initial  conditions  can  be 
specified  accurately  enough  for  meaningful  results  to  be  obtained 
from  the  calculation. 

It  has  been  demonstrated  that  the  choice  of  diffusion  model  can 
have  a  significant  influence  on  the  results  even  for  flows  with 
large  axial  vorticity.  This  is  in  agreement  with  an  identical 
comparison  for  2-D  mixing  calculations  with  zero  vorticity.  It  is 
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concluded  that  the  use  of  the  binary  diffusion  model  should  be 
abandoned  in  chemical  laser  mixing  calculations. 

Comparisons  with  the  premixed  model  indicate  that  the  use  of  the 
premixed  assumption  will  significantly  underpredict  the  amount  of 
iodine  dissociation  and  small  signal  gain.  This  result  should  put 
to  rest  the  misconception  that  the  premixed  assumption  always 
predicts  the  most  optimistic  laser  performance.  However,  it  has 
been  shown  that  the  use  of  the  premixed  model  with  the  use  of  an 
artificially  increased  I2  dissociation  rate  constant  will  give  the 
correct  average  values  for  yield  and  gain.  This  result  provides 
some  justification  for  the  use  of  the  premixed  model  in  situations 
where  the  degree  of  I2  dissociation  is  known  from  experiments. 

It  is  believed  that  the  TRIMIX  code  provides  a  rational  model  for 
analyzing  the  performance  of  oxygen/iodine  lasers  at  a  reasonable 
computational  cost.  It  includes  the  3-D  aspects  (geometry  and 
axial  vorticity)  which  are  ignored  in  2-D  codes  and  it  is  at  least 
one  and  possibly  two  orders  of  magnitude  faster  than  an  exact 
Navier-Stokes  code.  However,  the  code  does  contain  approximations 
in  the  definition  of  the  starting  conditions  and  it  should  be 
compared  Vvith  both  experimental  data  and  a  Navier-Stokes  solution. 
Future  effort  should  concentrate  on  adding  the  capability  to  solve 
the  species  equations  simultaneously,  addition  of  a  resonator  model 
for  power  extraction,  and  code  validation  efforts. 
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APPENDIX  A 

VORTICITY  TRANSPORT  EQUATION  FOR  A  COMPRESSIBLE , 
VARIABLE  VISCOSITY  FLUID 


The  vorticity  transport  equation  for  an  incompressible,  constant 
viscosity  fluid  has  a  particularly  simple  form  and  may  be  found  in 
many  books.  The  corresponding  result  for  nonconstant  viscosity  and 
a  compressible  fluid  is  difficult  to  locate  and  is  presented  in 
this  appendix. 

The  vorticity  and  velocity  are  defined  by: 


a  =  VxV  =  wxex  +  wyey  +  u)zez  (A.l) 

V  =  ue  +  ve  +  we 

A  Jf  * 

The  momentum  equation  for  steady  flow  may  be  written  as: 


-  pvxw  =  -  Vp  +  V  (pVV) 


(A. 2) 


By  taking  the  curl  of  Equation  A.  2  and  combining  with  the 
continuity  equation,  the  vorticity  transport  equation  is  obtained 
as: 


(pV-V)w 


(w-V)pV 


+  VpxV 


=  V(pVw)  -  (Vp-V)a/  +  Vx[(Vp*V)V]  +  VpxV  2 v 


(A. 3) 


Equation  A. 3  is  a  vector  relationship  which  provides  the  transport 
equations  for  the  scalar  components  wx,  wy  and  wz.  Note  that  for 
constant  density  and  viscosity,  the  gradients  of  p  and  p  vanish  and 
the  usual  equation  is  recovered. 
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Since  only  the  axial  vorticity  (ux)  is  solved  for  in  the  TRIMIX 
code,  we  express  Equation  A. 3  for  the  scalar  component  as 
follows: 


For  notational  convenience,  we  define: 

component  of:  VpxV  ^V2/2 


or  s  x 


(A. 4) 


P  =  x  component  of:  Vpxv  V  +  Vx[(Vp*V)V]  -  (Vp*V)w 


The  scalar  equation  for  ui  is  obtained  as: 


pV-Vw  +  a  =  V*  fpVw  +  wVpu  +  /? 

X  \  X) 


(A. 5) 


Equation  A. 5  is  parabolized  by  dropping  the  x  derivatives  in  the 
expression  for  P  and  the  x  derivative  in  the  divergence  of  pVwx. 

The  form  of  Equation  A. 5  which  is  solved  for  w  is  given  by: 


PV‘Vwx  =  V*  (pVwJ  +  (ox  f^(pu)  +  Q  (A.  6) 
Q  =  wy  fprCPu)  +  wz  f^(pu)  +  p  -  a 

The  vorticity  components  in  the  y  and  z  direction  in  the  expression 
for  Q  are  calculated  from  the  definitions  of  Equation  (A.l).  That 
is: 


du  <3w 
y  dz  ~  dx 


_  3v  du 
z  dx  ~  dy 


(A.l) 
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APPENDIX  B 
INITIAL  CONDITIONS 


In  this  appendix,  the  global  conservation  equations  for  mass, 
momentum  and  energy  are  used  in  conjunction  with  some  closure 
assumptions  to  determine  the  initial  flow  conditions  at  station  D 
(see  Figure  2.2). 

The  axial  momentum  equation  at  station  D  may  be  written  as: 


’pmA)D  -  Mp  +  (v)x  +  M2  f®-1) 

pm  s  p  +  P*2 

where  subscript  p  refers  to  the  primary  flow  and  subscript  1  and  2 
refer  to  the  injected  jets. 

The  flow  area  at  station  D  must  be  conserved  such  that: 


A  =  Ap  +  A^  +  A2  (B.2) 

Assuming  normal  injection,  such  that  the  jets  add  nothing  to  the 
axial  momentum,  and  ignorirg  wall  friction,  the  momentum  at  station 
0  and  D  is  equal.  That  it. 


( B .  3 ) 


Ignoring  heat  transfer  between  the  primary  flow  and  the  injected 
jets,  the  total  temperatures  at  station  D  are  known.  The  energy 
equation  for  the  jets  and  the  primary  are  given  by: 


1  +  (yp-i)Mp  /2 
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(B.4) 


« 


Tos  -  Ti  i1  +  (VS-D»12/2) 

Tos  *  T2  i1  +  (Vs-DM22/2) 

where  Top  is  the  primary  flow  total  temperature  at  station  0,  and 
Tos  is  the  plenum  temperature  for  the  jets.  Subscript  s  refers  to 
the  secondary  flow  (i.e.,  the  jets). 

Since  there  is  assumed  to  be  no  mixing  between  the  injection 
locations  and  station  D,  the  mass  flows  are  given  by: 


mp  =  (pUA)  p 

=  (pUA) x  (B. 5) 

m2  =  (pUA)  2 

Since  both  jets  are  fed  from  a  common  plenum,  the  jet  mass  flows 
are  related  by  the  ratio  of  the  injection  hole  areas.  That  is, 


(B.  6) 


With  the  assumption  that  the  flows  are  pressure  equilibrated  at 
station  D,  the  pressures  are  given  by: 


To  close  the  system  of  equations,  a  value  for  the  ratio  of  the  jet 
areas  at  station  D  must  be  assumed.  The  obvious  choice  is  to  use 
the  same  ratio  as  the  hole  pattern.  This  gives: 


VA2 


( B .  8  ) 
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The  given  information  consists  of  (PJOo'  Top,  Tos,  mp,  m,,  m2, 
(A.,/A2)  inj-  and  the  area  A  at  station  D.  For  given  values  of  A1/A2  and 
the  pressure  P,  all  the  flow  values  and  the  areas  occupied  by  the 
jets  and  primary  may  be  calculated  at  station  D.  As  indicated. 
Equation  (B.8)  is  used  for  the  jet  area  ratio.  The  pressure  P  is 
obtained  by  assuming  that  it  is  the  same  value  as  the  premixed 
pressure.  The  premixed  pressure  is  easily  calculated  by  assuming 
instantaneous  mixing  of  the  jets  and  primary  without  chemical 
reactions.  The  premixed  conservation  equations  for  mass,  momentum 
and  energy  are  easily  derived  and  are  not  repeated  here. 

To  solve  the  system  of  equations,  an  iterative  solution  is  required 
where  the  Mach  number  M1  is  the  iteration  parameter.  Note  that  the 
choice  for  the  jet  area  ratio  (Eq.  B.8)  dictates  that  M2  =  M1  and 
therefore  the  properties  (temperature,  velocity  and  density)  of 
both  jets  are  identical. 

To  summarize:  With  the  assumption  of  negligible  mixing  and  heat 
transfer  and  rapid  pressure  equilibration  between  the  jets  and  the 
primary,  the  global  conservation  equations  may  be  solved  to  obtain 
the  initial  condition  for  the  3-D  parabolic  mixing  calculation.  It 
is  believed  that  these  approximations  and  the  closure  assumptions 
of  premixed  pressure  and  Equation  (B.8)  for  the  jet  area  ratio  are 
reasonable.  To  assess  the  accuracy  of  these  approximations  would 
require  generating  some  exact  3-D  solutions  of  the  Navier-Stokes 
equations  for  comparison  purposes. 
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APPENDIX  C 
GAIN  MODEL 

The  gain  is  dependent  upon  the  density,  temperature  and  I  and  I* 
concentrations  and  is  obtained  from  Equation  C.20  and  C.17  of 
Reference  C.l. 

a  =  A2A21G  [N2  -  g^/g^/Srr  (C.l) 

where  A  is  the  wavelength  and  A21  is  the  Einstein  coefficient  for 
spontaneous  emission.  The  line  shape  function,  G,  is  obtained  from 
the  Voight  function  and  includes  both  the  effects  of  Doppler  and 
pressure  broadening.  Pressure  broadening  was  included  in  the  gain 
calculation  through  the  use  of  experimental  pressure  broadening 
coefficients  obtained  from  the  literature. 

Schlie  (Ref.  C.2)  reports  that  the  laser  should  operate  on  the  3-4 
hyperfine  transition  and  this  is  assumed  to  be  the  case  here.  For 
this  transition,  the  number  densities  of  the  upper  and  lower  states 
are  related  to  the  total  number  density  of  I  and  I*  by  the 
following  (Ref.  C.2): 

N,  =  9N,/24 

(C.2) 

N2  «  7Nj„/12 

where  the  number  density  of  any  species  is  related  to  its  mass 
fraction  by: 

Nj  =  Na  pKj/M,-  (C .  3 ) 

where  NA  is  Avogadro's  number  and  M,-  is  the  molecular  weight  of 
species  i. 

The  degeneracy  ratio  and  the  Einstein  coefficient  are  (Ref.  C.2): 
g2  =  7g,/9  A21  =  5.0  s'1  (C.4) 


57 


The  Voight  (Ref.  C.3)  line  shape  function  is  given  by: 

G  =  w~  ir^  [1  "  erf(y)]  EXp(y2) 

(C.5) 

y  =  Avl  /ln2/Av0 

where  AvL  and  Av0  are  the  Lorentz  and  Doppler  full-width,  half 
maximum  line  shape.  The  Doppler  width  is  given  by: 

Av0  =  2  ( 2RTln2/M)  1/2/A  (C.6) 

For  pressure  broadening,  the  line  width  (Eq.  4.8  of  Ref.  C.l)  is 
proportional  to  the  collision  frequency  (the  inverse  of  the 
average  time  between  collisions) .  Using  the  expression  of 
Reference  C.4  (page  1023)  for  the  collision  frequency,  the  Lorentz 
line  width  may  be  expressed  as: 


where  Tr  is  any  reference  temperature,  NA  is  Avogadro's  number,  and 
Ma  and  Mj  are  the  molecular  weights  of  species  a  and  i.  Qaj  is  the 
optical  cross  section,  and  P,.  is  the  partial  pressure  of  species  i. 

Generally,  an  expression  for  the  optical  cross  section  Qai  is  not 
available,  and  the  Lorentz  width  is  evaluated  from  measured 
pressure-broadening  coefficients.  The  literature  values  (denoted 
by  craj)  are  usually  presented  in  a  manner  which  implies  no 
temperature  dependence  on  AvL.  That  is, 

Avl  =  l  OfgjPj  (C.8) 
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The  /?aj  coefficients  are  related  to  the  measured  literature 
values  of  aai  by: 


where  Tm  is  the  measured  gas  temperature  at  which  oraj  was 
determined. 


If  the  optical  cross  section  is  independent  of  temperature,  then  n 
=  0.50  from  Equation  C.7.  On  page  1025  of  Reference  C.4  it  is 
stated  that  Qaj  is  inversely  proportional  to  the  square  root  of 
temperature.  If  so,  then  it  follows  that  n  =  1.0.  If  the  measured 
ocai  are  available  at  room  temperature,  then  by  choosing  the 
reference  temperature  equal  to  room  temperature  (i.e.,  Tr  =  295  K) 
we  have  (3ai  =  araj .  This  is  the  approach  used  in  the  code :  the  /?ai 
are  taken  to  be  the  measured  room  temperature  values  and  Equation 
C . 7  is  used  to  calculate  the  Lorentz  width. 
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